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We calculate the change in susceptibility resulting from a thin sheet with reduced 
penetration depth embedded perpendicular to the surface of an isotropic supercon- 
ductor, in a geometry applicable to scanning Superconducting QUantum Interference 
Device (SQUID) microscopy, by numerically solving Maxwell's and London's equa- 
tions using the finite element method. The predicted stripes in susceptibility agree 
well in shape with the observations of Kalisky et a/.[r of enhanced susceptibility 
above twin planes in the underdoped pnictide superconductor Ba(Fei_2;Coa;)2As2 
(Ba-122). By comparing the predicted stripe amplitudes with experiment and using 
the London relation between penetration depth and superfluid density, we estimate 
the enhanced Cooper pair density on the twin planes, and the barrier force for a 
vortex to cross a twin plane. Fits to the observed temperature dependence of the 
stripe amplitude suggest that the twin planes have a higher critical temperature than 
the bulk, although stripes are not observed above the bulk critical temperature. 



I. INTRODUCTION 



Scanning SQUID microscopy has been used very successfully to elucidate fundamental 
properties of superconductors from spatially resolved magnetometry images, including the 
determination of the pairing symmetry of the cuprate superconductors [2j , tests of the inter- 
layer tunneling model for high temperature superconductivity [3J, and the placement of limits 



2 



on spin-charge separation in high temperature superconductors [1] . Recently a new dimen- 
sion has been added to scanning SQUID measurements: scanning SQUID susceptometry [H] 
has enabled spatially resolved measurements of superconducting penetration depths [SI E] , 
the observation of spontaneous persistent currents in single mesoscopic normal rings [8j , and 
measurements of fluctuations in single mesoscopic superconducting rings [9]. An advantage 
of scanning SQUID susceptibility measurements, aside from their exceptional sensitivity, is 
that they can be easily, reliably, and quantitatively modeled [10]. Recently Kalisky et al. [1] 
reported stripes of enhanced susceptibility in scanning SQUID susceptometer measurements 
of single crystals of the pnictide superconductor Ba(Fei_a;Co2.)2,.As2. They associated these 
stripes with twin boundaries in the superconductor, for a number of reasons: 1) The stripe 
spacings and orientation were consistent with optical and x-ray observations of stripes in 
the same ^ and similar [11] samples. 2) The stripes were only observed for underdoped 
samples, which undergo a tetragonal to orthorhombic crystal structure transition and there- 
fore have twins, but not for optimally doped or overdoped samples. 3) The stripes changed 
position when the samples were warmed above the tetragonal-orthorhombic transition tem- 
perature, but not when the sample was only warmed above the superconducting transition 
temperature. Further, it was observed that vortices did not pin on the stripes, and that 
when dragged using either a SQUID sensor, or a magnetic force microscope tip, the vortices 
did not cross the stripes. Enhanced diamagnetic susceptibility and enhanced critical tem- 
peratures associated with twinning have been reported previously from bulk measurements 
of several elemental superconductors [T^ . 

The observation of stripes in susceptibility is quite interesting qualitatively, because it 
indicates that the superconductivity is different on the twin boundaries than in the bulk 
in these novel superconductors. However, in order to fully understand these results it is 
important to model them quantitatively. This is a challenge, because of the special sample 
geometry involved. Kogan [TO] developed a theory for the Meissner response of anisotropic 
superconductors to several types of locally applied magnetic fields, including from a circular 
field coil. This theory produces exact solutions for the problem of scanning SQUID suscep- 
tometry of a homogeneous superconductor. However, it is not immediately apparent how 
to apply this theory to our geometry. In the present paper we use numerical methods to 
solve the problem of local susceptibility measurements of an inhomogeneous superconductor 
with a planar sheet with different superconducting properties oriented perpendicular to the 
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FIG. 1: Scanning SQUID susceptibility image of an underdoped (j;=0.051, Tc=18.25K) single 
crystal of Ba(Fei_2,'Cox02As2 at T=17K. The arrow indicates the direction along which the stripes 
were averaged to form the cross-section plotted in Figure [2] 

bulk surface. We apply the results of this modeling to the experiments of Kalisky et al. 
and find good agreement. Our analysis indicates that there is substantial additional Cooper 
pair density on the twin boundaries, and that the critical temperature of the twin boundary 
region is higher than that in the bulk. 

II. ENHANCED SUSCEPTIBILITY IMAGES 

An example of the stripes of enhanced susceptibility observed by Kalisky et al. is shown 
in Figure [1} A scanning SQUID susceptometer sensor has a single turn field coil surrounding 
a co-planar pickup loop integrated into the SQUID sensor. Susceptometry measurements 
are made by applying a small alternating current to the field coil and phase sensitively 
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FIG. 2: (a) Averaged cross section along the stripes direction of the bottom 3/5 of the data of 
Figure [T] (symbols, dashed line). Subtracting out the piece- wise linear background (dotted line in 
(a)) results in the symbols and dashed line in (b). The solid line in (b) is modeling as described in 
the text. 

sensing the response of the SQUID, proportional to the flux through the fleld coil, to this 
current. The data shown in Figure [l] was taken at T=17 K using a sensor [13] with an 
effective fleld coil radius of R=8.85fim and an effective pickup loop radius of r=2.1/im. The 
effective height of the sensor above the sample surface was ZQ=1.5fim, derived from fltting 
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magnetometry images of individual superconducting vortices in the sample. The ac current 
through the SQUID was 0.25 mA, corresponding to a magnetic induction at the sample 
surface of 17/iT. The sign convention is chosen such that diamagnetic shielding is positive: 
higher numbers and white colors represent stronger diamagnetic shielding. 

In order to compare experimental data with our modeling, we averaged the image in 
Figure [T] along the direction indicated by the arrow to obtain the cross-section displayed 
in Figure [2^. There are broad spatial variations in the susceptibility in addition to the 
stripes. We subtract them from the data using a piece-wise linear background indicated 
by the dashed line in Fig. [2^. This results in the data displayed as symbols in Fig. |2}d. 
The susceptibility peaks in this cross-sectional average have amplitudes of 3.37 ± 0.35$o/^ 
and full widths at half-maximum of 8.92 ± 0.62yum. This is to be compared with a low 
temperature susceptibility of about 600 ^o/A. The solid line in Figure [2] is the result of 
modehng as described below. 

III. NUMERICAL MODEL 

We believe that the stripes in susceptibility seen in Fig. 's [T] and [2] are due to reduced 
penetration depths associated with the twins. We have not been able to perform an ana- 
lytical calculation of the change in susceptibility due to a narrow sheet of superconductor 
with reduced penetration depth. Instead, we numerically solved the coupled Maxwell's and 
London's equations in the problem with COMSOL, a commercial program that uses finite 
element methods to solve partial differential equations. The geometry we assumed is shown 
in Figure [3] We used COMSOL to numerically solve Laplace's equation V^i? = for 
z > 0, and Londons' equation V^if = H/\^ for z < 0, with A = in the sheet volume 
\x — Xq\ < w /2, 2; < 0, and A = A^ in the volumes \x — Xq\ > w/2, ^ < 0. The boundary con- 
ditions used were continuity of H across the plane z = and the half-planes x = xq — w/2, 
and X = xq + w/2 {z < 0), H = on the outside surfaces of the enclosing space, and 
Ht = I /2TTa, where Ht is the component of H tangential to the surface of the torus along 
the direction of maximum surface curvature. This last boundary condition was derived by 
using Ampere's law along circles around the torus walls. In general the magnetic fields 
are not constant along these circles, but this approximation becomes exact as a — 0. For 
these simulations we used a/R = 0.05. Doubling a changed the computed results by only 
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FIG. 3: Modeling geometry: (a) 3-d rendering; (b) top view; (c) side view; and (d) expanded view 
of the field coil. A superconductor occupying the half-space z<0 is assumed to have a penetration 
depth Xf, except for a sheet of width w with penetration depth As centered at x = xq- The 
susceptometer field coil, represented by a torus with major radius R and minor radius a oriented 
with its plane parallel to the surface of the superconductor and centered at x = 0,y = 0,z = zq, 
carries a current /. The susceptibility x = ^/^ calculated by integrating the z-component of 
the resultant magnetic field over a circle of radius r centered at (0, 0, zq) to obtain the flux ^. For 
the modeling in this paper we used -R=8.85^m, zq/R=0.17, a/i?=0.05, r/i?=0.25. The volume of 
space modeled was of size 12x12x8 R^. 

a few percent. The simulations shown here used meshes generated by an advancing front 
algorithm with the COMSOL setting hauto, which automatically sets several parameters 
for the meshing, equal to 4. This resulted in about 600000 degrees of freedom. Using the 
COMSOL mesh density setting hauto = 4, with 4 times fewer degrees of freedom, resulted 
in values for the 2;-component of the magnetic field at (0, 0, zo) which were different from 
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FIG. 4: Calculated normalized supercurrent densities at z=0 for w/R=0.2 in the x (a,c) and y 
(b,d) directions for the sheet penetration depth Xs/R=0.2 and bulk penetration depth Afe/i?=0.2. 
(a,b), and As/i?=0.1, Xb/R=0.2 (c,d). The small scale inhomogeneities in the image are due to 
discretization error. 

those for the more dense mesh by a few percent. For all the modeling presented here we 
used the denser mesh, which took about 8 minutes per point on a 2.5 GHz Intel Core 2 Duo 
Mac Book Pro. 

IV. RESULTS 

Figure [2] compares the experimental cross-section from Fig. [T] with our numerical re- 
sults, obtained by integrating H^, the calculated field in the z direction, using w/R = 0.2, 
Xs/R=0.1, and Aft/-R=0.12, over the area of the pickup loop for a number of different val- 
ues of Xq. Noise is visible in the modeling due to the small difference between and 
required to fit the experimental stripe heights, given the value for w/R chosen. Smaller 
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values of w/R require more mesh elements than can be accommodated by our computer 
memory. We will discuss how the stripe peak heights scale with the different parameters 
below. The susceptibility contributions from the 8 stripes visible in the image were added. 
The fitting parameters in this plot were Xb, Xs, w, the center positions of the stripes and 
a uniform vertical shift. This figure shows that the stripe peak shapes and widths agree 
with experiment within the variation from stripe to stripe: The simulated peak heights are 
3.43 ± 0.01$o/^ while the experimental peak heights are 3.37 ± 0.35$o/^- The experi- 
mental peak full-widths at half-maximum (FWHM) are 8.92 ± 0.62/im, while the simulated 
peak widths are 9.34 ± 0.82/im. The stripe susceptibility peak widths are limited by ex- 
perimental resolution rather than their intrinsic widths: Equally good agreement with the 
experimental peak width can be obtained for any value of w below about 5/xm. Increasing 
the simulated width to w/R = 0.6, corresponding to w = b.Sfim, gives a simulated FWHM 
of 10.07 ± 0.75/im, about one standard deviation larger than the experimental peak widths. 

Figure |4] shows results for the dimensionless current densities j^R'^/I and jyR'^/I (using 
j = V X i?) in the x, y directions respectively at z = 0, the surface of the superconductor, 
induced by the field coil. In the case in which the penetration depth of the sheet is the same as 
the rest of the superconductor (Fig. |4^,b, Xs/R = XbR=0.2, ) the field coil induces screening 
currents which are strongest directly under the ring, and circulate with the opposite sense 
as the currents in the field coil. If the sheet penetration depth is smaller than that of the 
bulk (Fig. |4] c,d, Xs/R=0.1, Xb/R=0.2) there is an additional component of the screening 
current in the x direction concentrated under the field coil, and a more delocalized excess 
component in the y direction. We chose a larger difference between A;, and for this image 
than for the fit of Fig. [2] to increase the contrast in the sheet region. 

V. SCALING 

As discussed above, the experimental width of the stripes is resolution limited. The 
sheets of enhanced superfiuid density could be as narrow as the coherence length (~3-4 
nm [Hj). We therefore performed simulations to see how the predicted results scaled as w 
became small. The results are shown in Figures |5]j6j In Figure [5] we plot the dimensionless 
susceptibility peak height AHzR/I, the difference between the ^-component of the magnetic 
field at the center of the field coil induced by the current I for xq = minus that at xq = 2R, 



9 



0.06 
0.05 



DC 

I 

I 

< 



0.04 
N 0.03 
0.02 
0.01 











/ / 



+ X O 



(b) 0.1! 



0.1 



DC 

^ 

< 



0.05 



10 



0. 



o 



/ + 







0.5 

(w/R) 



1/2 



FIG. 5: Dependence of the stripe peak height l^HzR/ 1 vs. 1/As(a) and lU"'^/^ (b), where Aff^ is the 
difference between xq / R = and xq/R = 2 oi the z-component of the magnetic field at the center 
of the field coil induced by the field coil current /. The symbols are the numerical simulation and 
the dashed lines are linear fits. 

as a function of R/\s (a) and vs. (w/RY^'^ (b). The dashed hnes in this figure are hnear 
fits. The susceptibihty peak heights are roughly proportional to l/A^, with a slope which is 
proportional to A;, — A^ (a), and roughly proportional to w^^'^ (b). The two figures Fig. [s^ 
and jsj) indicate that the susceptibility stripe amplitudes scale as w^^'^/Xg, proportional to 
the square root of the two-dimension sheet Cooper pair density. In Figure |6] we have plotted 
the dimensionless normalized susceptibility peak height XsHzR^/"^ / Iw^/"^ versus (Ab — Xs)/R, 
the reduced difference in penetration depths between the stripe and the bulk. The scaling 
works reasonably well over the range of parameters chosen. The solid line in Figure [6] is the 
empirical relation 

\ AH Ri/2 / X, _ X \ 

(1) 



XsAH.R^/^ , f^Xb-Xs 

= a tanh p- 



with a = 0.021 and f3 = 2.3796. We use this relation in the modeling that follows. 
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FIG. 6: Plot of the dimensionless stripe susceptibility peak height XsAHzR^/'^ / Iw^^"^ vs {Xfy — Xg) / R 
for various values of the reduced sheet width w/R and sheet penetration depth Xg/R. For large 
values of Af, — As the susceptibility peak height is nearly independent of Xf,, indicating that the 
susceptibility peak height is proportional to w^^'^/Xg, and therefore proportional to the square root 
of the two-dimensional sheet Cooper pair density Ng. The susceptibility peak height approaches 
zero linearly as Xh Xg. The solid line is the empirical function Eq. [T] used for the modeling of 
the temperature dependence of the stripe amplitude below. 



VI. TEMPERATURE DEPENDENCE 



Kalisky et al. found that as the sample temperature approached Tc the stripe peak 
amplitudes became larger. Experimental data on an underdoped (x=5.1%) sample is shown 
as the symbols in Figure [7} The solid symbols in this figure were obtained by fitting an 
image with multiple stripes to our numerical model, using the procedure of Fig. [2| but with 
a single vertical scaling factor, for a number of sample temperatures. The vertical error bars 
were obtained by varying the scaling factor from its optimal value until the between fit 
and experiment doubled. We estimate that the errors in temperature are smaller than the 
symbol widths. Our modeling shows that if the sheet and bulk penetration depths had the 
same temperature dependence there would be no change in the stripe amplitude until the 
bulk penetration depth became comparable to the radius of the field coil - very close to T^. 
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FIG. 7: Plot of the temperature dependence of the stripe amphtude Ax(T) (symbols). The solid 
lines are best fits to the data using the empirical formula Eq. [T] and the two-fluid temperature de- 
pendence (Eq. [2j p=4) for the penetration depths, and different critical temperatures T^s and T^^j, 
for the sheet and bulk respectively, for widths w/R = 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, cor- 
responding to the sheet width w between 4.4/im and 17.6nm. The best fit values and uncertainties 
for As and Tc^s/T^fi are plotted in Figure |8j 

We have considered two scenarios for a difference in temperature dependences in the stripe 
amplitudes: 1) The bulk and sheet critical temperatures are the same, but the penetration 
depths have different temperature dependences below T^. Such a difference could occur, for 
example, if the bulk is a weak coupling superconductor, and the sheet is a strong coupling 
one, or visa versa. In addition, such a difference could occur if the pairing symmetry is 
different in the sheet than in the bulk. The temperature dependence of the penetration 
depth can be parameterized as [15] 

A(T) = A(0)/v'l - {T/T,y (2) 

where, for example, good fits to the BCS predictions can be obtained using p=2 for s-wave 
and p=4/3 for (i-wave pairing symmetries [l5j. We fit the experimental data of Fig. [T] using 
the Eq. [Tjand Eq. |2]using A^, Ps, and ph (the p exponents for the sheet and bulk respectively) 
as the three fitting parameters. In all cases the fits gave unreasonably large values for the 
p's. For example, the best fit assuming w/R=0.2, Ab=0.04 gave Ps=34.2 and Ph=l7.9. It 
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FIG. 8: Best fit values and uncertainties (using a doubling of the best fit value for as a criterion) 
for the parameters \s/R and Tc^s/Tc^i) obtained by fitting the data of Figure [Tj for various assumed 
values for the sheet width assuming \b{T = 0)/i? = 0.04. 

therefore seems unlikely that the first scenario, in which the twin boundaries and bulk have 
the same critical temperature, can explain the data. 

2) On the other hand, good fits to the data can be obtained by assuming different critical 
temperatures for the bulk and sheet. The solid lines in Figure [7] are a superposition of fits 
of Eq. [l] to the data, assuming that both the sheet and bulk penetration depths follow the 
two-fluid temperature dependence {p = 4) in Eq. [2| assuming two different Tc's for the 
bulk and the sheet, for various values of w/R, and Af,(T = 0)/-R=0.04 (corresponding to 
Afe = 0.35/im[TB]). The quality of the fits is not significantly different for different values 
of the sheet width w within the physically reasonable range. Similar results, with slightly 
different values for Tc^s/Tc^, are obtained using p = 2 (BCS) ot p = 4/3 (rf-wave) in Eq. 
|2] for the fits. Figure [s] displays the best fit values and error bars (using a doubling of the 
best-fit as the criterion) for Xs/R and T^^s/Tcfi from these fits. In this modeling we assume 
that w is independent of temperature. Since the stripe amplitudes scale like w^^^, any w 
temperature dependence can be neglected relative to the temperature dependence of A;, — A^. 
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FIG. 9: Best-fit values, and upper and lower bounds for the enhancement of the two-dimensional 
sheet Cooper pair density as a function of sheet width from fits to experimental susceptibility stripe 
amplitudes. 

VII. DISCUSSION 

Our analysis above indicates that the twin boundaries in underdoped Ba(Fei_2,Coa;)xAs2 
have a shorter penetration depth than the bulk. A shorter penetration depth is usually 
associated with a higher superfiuid density, although the relation between the two quantities 
could be modified by fluctuations [17] . We assume for simplicity the London relation to 
obtain an estimate for the excess 2-dimensional sheet Cooper pair density ANg 

m* / 1 1 \ 

ANs = W——r - T2 (3) 

where m* is the Cooper pair mass (assumed to be twice the bare mass of the electron) and 
e*=2e is the Cooper pair charge. From Fig. [9] we estimate that the excess 2-dimensional 
sheet Cooper pair density is lO^^m^^ < AN^ < 10^°m~^. 

Our analysis also indicates that the sheets have a higher critical temperature than the 
bulk, although stripes have not yet been observed above the bulk critical temperature, 
possibly because of superconducting fluctuations or local variations in the bulk Tc. 

Kalisky et al. also found that vortices were not pinned on the sheets, and that it was 
not possible to drag vortices across the sheets. The energy required to form a vortex in a 
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superconductor is given approximately by [18]: 



E--JZ^2^^'^^ (4) 



47r/ioA^ 

where k ^ A/^ ~ 140 [19], and L is the vortex length. From the estimates of Figure [s] 
we find 1 X lO^^^J < < 4 x lO^^^J for the difference in energy of the vortex on vs. 
off the sheet, assuming the crystal thickness L = lO^m. These energies are much larger 
than ksTc. These are very rough estimates, because it may well be that the sheet width is 
much smaller than the penetration depth, in which case much of the vortex field energy is 
outside of the sheet. We estimate from dragging experiments with an MFM tip that a force 
of 40 pN at 5K and 6pN at 14K was not enough to make a vortex cross a sheet. Using w 
as a characteristic length, the maximum and minimum excess vortex energies on the sheet 
correspond to forces Fg ~ AE^/w of 1.4 x 10~^A^ > > 9.2 x 10^^^ A^, easily large enough 
to explain the inability to drag vortices across the twin boundaries. 

Subsequent to the work of Kalisky et al. [1], Prozorov et al. [20j noted an enhancement 
of the critical current of slightly underdoped single crystals of Ba(Fei_^Co2;)2As2, which 
they associated with twin boundaries. The measurements of Kalisky et al. provide two 
mechanisms for this critical current enhancement: 1) Enhanced superfiuid density along the 
twin boundaries provides channels with enhanced depairing current densities, and 2) The 
enhanced superfiuid densities in the twin boundaries provides barriers to transverse vortex 
motion. Our modeling shows that the latter effect can be substantial. 



VIII. CONCLUSIONS 



In conclusion, we have used finite element methods to numerically solve Maxwell's and 
London's equations for the problem of enhanced susceptibility from a narrow sheet of super- 
conductor with reduced penetration depth embedded in a bulk superconductor in a geome- 
try appropriate for scanning SQUID susceptometry measurements. We find good agreement 
between our modeling and experiment for the lineshape for cross-sections across stripes in 
enhanced susceptibility measurements on underdoped samples of the pnictide superconduc- 
tor Ba-122. By scaling our simulations and comparing the results with experiment we obtain 
estimates of the enhanced Cooper pair sheet density on the sheet. The barrier energies and 
forces we estimate are large enough to explain the observed lack of pinning on the sheets. 
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and the experimental inability to drag vortices across the sheets. 
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